#This file creates main findings for Dahlgaard (2016)

load("het_dat_last.rdata")

## commented out for replication files 

#  install.packages("xtable")
#  install.packages("ggplot2")
#  install.packages("rmeta")
#  install.packages("dplyr")
#  install.packages("rdrobust")

library(xtable)
library(ggplot2)
library(dplyr)
library(rmeta)
library(rdrobust)

## Estimate  models using data split on year and covariates
mean_eff <- array(NA, dim = c(4, 2) )
mean_se  <- array(NA, dim = c(4, 2) )
mean_eff_prob <- array(NA, dim = c(4, 2) )
upp_prob  <- array(NA, dim = c(4, 2) )
low_prob  <- array(NA, dim = c(4, 2) )
years <- c(2013, 2014)

rownames(mean_eff) <- c("2009-2013",
                        "2013-2014",
                        "2009-2013 rob",
                        "2013-2014 rob")
colnames(mean_eff) <- c("Abstained", "Voted")

rownames(mean_se)  <- rownames(mean_eff)
colnames(mean_se)  <- colnames(mean_eff)

rownames(mean_eff_prob)  <- rownames(mean_eff)
colnames(mean_eff_prob)  <- colnames(mean_eff)
rownames(upp_prob)  <- rownames(mean_eff)
colnames(upp_prob)  <- colnames(mean_eff)
rownames(low_prob)  <- rownames(mean_eff)
colnames(low_prob)  <- colnames(mean_eff)

for (k in 1:2){
  data <- 
    data.frame(data_list_votes[[1]]) %>%
    filter(voted_last == k - 1) %>%
    mutate(treated = two.days > 0)
  
  
  data_close <- 
    data %>%
    ungroup() %>%
    mutate(voted     = round(par_vote*n),
           abstained = round((1-par_vote)*n)) %>% 
    select(year, two.days, treated, voted, abstained)
  
  data_voted <- 
    as.data.frame(lapply(data_close, 
                         function(x,p) rep(x,p), 
                         data_close[["voted"]])) %>%
    mutate(voted = 1)
  
  data_abstained <- 
    as.data.frame(lapply(data_close, 
                         function(x,p) rep(x,p), 
                         data_close[["abstained"]])) %>%
    mutate(voted = 0)
  
  data_master <- 
    rbind(data_voted, data_abstained)
  
  for(j in 1:2){
    sub_data <- 
      data_master %>%
      filter(year == years[j])
    mod <- 
      with(sub_data, rdrobust(y = voted, x = two.days))
    
    mean_vote <- with(subset(sub_data, two.days < 0), mean(voted))

        mean_eff[j,k] <- mod$coef[1]
    mean_se[j,k]  <- mod$se[1]
    mean_eff_prob[j,k] <- 
      qnorm(mod$coef[1] + mean_vote) - 
      qnorm(mean_vote)
    upp_prob[j,k] <- 
      qnorm(mod$coef[1] + mean_vote + qnorm(0.975)*mod$se[1]) - 
      qnorm(mean_vote)
    low_prob[j,k] <- 
      qnorm(mod$coef[1] + mean_vote + qnorm(0.025)*mod$se[1]) - 
      qnorm(mean_vote)  
    
    mean_eff[j + 2,k] <- mod$coef[3]
    mean_se[j + 2,k]  <- mod$se[3]
    mean_eff_prob[j + 2,k] <- 
      qnorm(mod$coef[3] + mean_vote) - 
      qnorm(mean_vote)
    upp_prob[j + 2,k] <- 
      qnorm(mod$coef[3] + mean_vote + qnorm(0.975)*mod$se[3]) - 
      qnorm(mean_vote)
    low_prob[j + 2,k] <- 
      qnorm(mod$coef[3] + mean_vote + qnorm(0.025)*mod$se[3]) - 
      qnorm(mean_vote)  
    }
  

}


# Plotting ----------------------------------------------------------------

mean_eff <- mean_eff * 100
mean_se  <- mean_se  * 100

effs <- mean_eff[1,]
low  <- mean_eff[1,] + mean_se[1,] * qnorm(0.025) 
upp  <- mean_eff[1,] + mean_se[1,] * qnorm(0.975) 


for (i in 2:4){
  effs <- c(effs, mean_eff[i,])
  low  <- c(low, mean_eff[i,] + mean_se[i,] * qnorm(0.025) ) 
  upp  <- c(upp, mean_eff[i,] + mean_se[i,] * qnorm(0.975) )
  
}


for (i in 1:4){
  effs <- c(effs, mean_eff_prob[i,])
  low  <- c(low , low_prob[i,])
  upp  <- c(upp , upp_prob[i,])
}


labels <- c("Effect in 2013 conditioned
            \non voting in 2009",
            "Effect in 2014 conditioned
            \non voting in 2013")

place <- rep(c(4.5, 3.5, 2, 1), 4)
voted <- rep(c("Abstained", "Voted"), 8)
type  <- rep(rep(c("Conventional RD", "Robust RD"), each = 4), 2)
prob  <- rep(c("Effect size in percentage points", "Effect size in probits"), each = 8)


plot_data <- 
  data_frame(effs, upp, low, place, voted, type, prob)

plot_prob <- 
  ggplot(data = plot_data,
         aes(y = place,
             x = effs,
             xmin = low,
             xmax = upp,
             fill = voted)) +
  geom_errorbarh(aes(alpha = voted, linetype = voted), height = 0) +
  facet_grid(type ~ prob, scales = "free_x") +
  scale_y_continuous(breaks = c(4, 1.5), labels = labels) +
  theme_classic() + 
  theme(axis.text = element_text(size = 12),
        axis.title = element_text(size = 20),
        legend.text = element_text(size = 16),
        legend.title= element_blank(), 
        legend.position = "top",
        panel.spacing = unit(2, "lines")) +
  ylab("") + 
  xlab("") +
  geom_point(aes(alpha = voted), size = 3) +
  scale_alpha_discrete(range = c(0.5, 1)) +
  geom_vline(xintercept = 0, linetype = 3, alpha = 0.5)

# table -------------------------------------------------------------------

tab <- matrix(NA, nrow = 4, ncol = 2)

for (i in 1:2){
  for (k in 1:2){
    data <- 
      data_list_votes[[1]] %>%      
      filter(voted_last == k - 1 & two.days < 0)
    data <- data %>% filter(year == years[i])
    tab[2*i - 1, k] <- paste(round(weighted.mean(data$par_vote, data$n)*100, 2)) 
    tab[2*i    , k] <- paste("(",sum(data$n),")", sep = "")
  }
}

colnames(tab) <- c("Abstained last", "Voted last")
rownames(tab) <- c("2013 conditioned on 2009",
                   "1",
                   "2014 conditioned on 2013",
                   "2")

xtable(tab)
